Ctrl+K

Методы стохастической оптимизации. Настройка гиперпараметров

Методы стохастической оптимизации. Настройка гиперпараметров#

from sys import platform
import sys
from matplotlib import pyplot
pyplot.rcParams['figure.dpi'] = 150
pyplot.rcParams['savefig.dpi'] = 150
if platform == "linux" or platform == "linux2" or platform == "darwin":
    sys.path.append("../../")

import random
import time
from matplotlib import pyplot

pyplot.rcParams['figure.dpi'] = 200
pyplot.rcParams['savefig.dpi'] = 200
from sys import platform
import sys

if platform == "linux" or platform == "linux2" or platform == "darwin":
    sys.path.append("../../")
%matplotlib widget
import numpy as np
from src.v2.impl.conditions import StepCountCondition, PrecisionCondition, AbsolutePrecisionCondition
from src.v2.impl.methods import CoordinateDescent
from src.v2.impl.metrics import StepCount, CallCount, GradientCallCount, HessianCallCount, ResultValue, RAMSize, ExecutionTime
from src.v2.impl.oraculs import LambdaOracul
from src.v2.runner.runner import Runner, FULL_VISUALIZE, VISUALIZE
from src.v2.visualization.animation import Animator
from src.v2.runner.runner import TABLE
from src.v2.impl.methods import GradientDescent, ScipyMethod, SimulatedAnnealing, Evolution, EvolutionRecombination, EvolutionCMA
from src.v2.impl.schedules import step_schedule, pow_step_schedule
from IPython.display import display, HTML
from math import exp, sqrt, cos, sin, pi, e

display(HTML("<style>pre { white-space: pre !important; }</style>"))


def seed(val: int = -1):
    if val == -1:
        val = time.time_ns()
    random.seed(val)
    np.random.seed(val % 2 ** 32)
    print("SEED:", val)


modules = [ExecutionTime(), RAMSize(), Animator(), StepCount(), CallCount(), GradientCallCount(), HessianCallCount(),
           ResultValue()]

Постановка задачи#

Главной задачей нашей работы является изучение теоретической и практической частей стахастических методов оптимизации. Для каждого метода будет представлено теоретическое обоснавание работоспособности метода. Чтобы определить эффективность методов в практическом применении будут рассмотрены задачи из предыдущих лабораторных работ, а также будет проведено сравнение различных метрик как и между стахастическими методами, так и с методами из предыдущих лабораторных работ.

Основное задание#

Метог отжига#

Первым методом, который будет предложен для изучения - метод имитации отжига. Как понятно из названия, данный метод основывается на имитации физического процесса. Суть метода заключается в случайной генерации следующего состояния системы на основе предыдущего. Если состояние оказывается лучше текущего, то совершается переход, если хуже, то переход совершается с каким-то шансом, зависящим от номера итерации. Обычно шанс на переход в худшее состояние имеет линейную убывающую зависимость, сами параметры очень сильно влияют на поведение метода, поэтому подбираются под каждую задачу отдельно.

Главное преимущество такого алгоритма - возможность выбираться из локальных минимумов, так как оказавшись в таком есть шанс несколько раз подряд оказаться в худшем состоянии и оказаться в области, откуда можно добраться до глобального минимума или другого локального. Сначала можно подумать, что тогда метод будет вести себя очень случайно и даже оказавшись возле глобального минимума будет иметь возможность выйти оттуда. Но если вспомнить что шанс на переход в худшее состояние - убывающая функция, то предположив, что метод постепенно приближался к глобальному минимуму, то уменьшив свой шанс на переход в худшее состояние до необходимого значения, он перекроет себе возможность на возвращение в те области из которых он пришёл и чьё состояние хуже текущей области.

В нашей реализации метод параметризуется с помощью learning rate, максимального отклонения от предыдущего состояния по каждой из осей, и decay, значения на которое уменьшается ганс на переход в худщее состояние после каждой итерации.

Чтобы сравнить метод отжига с другими методами рассмотрим авторскую двумерную функцию с большим количеством локальных минимумов. В качестве методов для сравнения будут выступать покоординатный спуск и Ньютоновский метод из библиотеки scipy.

seed(1717697282242590600)

methods = [
    SimulatedAnnealing(3, 0.7, 0.01),
    CoordinateDescent(learning_rate=3, aprox_dec=1e-30),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(100),
    AbsolutePrecisionCondition(1e-5, np.array([0]))
]

hard_function = LambdaOracul(
    lambda x: -5 * exp(cos(x) ** 2 - abs(x) / 6) + abs(x) / 5 * (sin(1.5 * x) * cos(2 * x / 3)) ** 2 + abs(sin(10 * x)))

result = Runner.run(methods, [hard_function], np.array([9.0]), conditions + modules, precision=1e-7,
                    **TABLE, **VISUALIZE, animate=True, animate_main=True, animation_main_full=True,
                    animator_main_step=0.01)

SEED: 1717697282242590600

LambdaOracul
+------------------------+-----------------+------------+-------------+-------------+---------------------+--------------------+------------+
| Method name            |   ExecutionTime | RAMSize    |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount |     Result |
+========================+=================+============+=============+=============+=====================+====================+============+
| Annealing((3, 0.7))    |      0.00194001 | 0.15625 Mb |          99 |         198 |                   0 |                  0 | -0.0338874 |
+------------------------+-----------------+------------+-------------+-------------+---------------------+--------------------+------------+
| CoordinateDescent(3)   |      0.00196099 | 0.0625 Mb  |          99 |         240 |                   0 |                  0 | 11.6239    |
+------------------------+-----------------+------------+-------------+-------------+---------------------+--------------------+------------+
| ScipyMethod(newton-cg) |      0.00746107 | 0.84375 Mb |          99 |          23 |                  11 |                  2 | 10.0645    |
+------------------------+-----------------+------------+-------------+-------------+---------------------+--------------------+------------+

В данном тесте метрики занимаемой памяти и времени исполнения немного подвели, так как и имитация отжига и координатный спуск отрабатывают так быстро, что сборщик мусора не успевает очищать память. Но не смотря на это видно превосходство метода имитации отжига - это единственный метод, который сошелся. Но из анализа результатов и наблюдений во время проведения тестов можно сделать вывод о невозможности достижения высокой точности этим методом в чистом виде. Так как область в которой случайным образом выбирается следующее состояние статична, то оказавшись возле минимума метод может лишь надеяться на удачу и пытаться оказаться ближе к минимуму. Помочь с этой проблемой может schedule, то есть понижание learning rate с увеличением числа итераций.

Чтобы определить эффективность добавления schedule для метода имитации отжига, предлагается провести сравнение на той же функции. В качестве schedule функции будет выступать алгоритм уменьшения learning rate в зависимости от количества числа итераций из предыдущей лабораторной работы, но с модификацией - теперь число итераций может возводиться в произвольную степень. Такая функция позволяет выставив очень маленький коэффициент и неединичную степень и таким образом почти не воздействовать на learning rate на первых итерациях. Например в нашем случае выбраны коэффициент 1e-6 и степень 3, значит на 100 шагу learning rate уменьшится всего в 2 раза, а на 500 шагу в 126 раз.

seed(1717700296282255800)

methods = [
    SimulatedAnnealing(3, 0.7, 0.01),
    pow_step_schedule(3.5, 0.000001, 3, SimulatedAnnealing)(3.5, 0.7, 0.01)
]

conditions = [
    StepCountCondition(500)
]

result = Runner.run(methods, [hard_function], np.array([9.0]), conditions + modules, precision=1e-7,
                    **TABLE, **VISUALIZE, animate=True, animate_main=True, animation_main_full=True,
                    animator_main_step=0.01)

SEED: 1717700296282255800

LambdaOracul
+--------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| Method name                    |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount |       Result |
+================================+=================+=============+=============+=============+=====================+====================+==============+
| Annealing((3, 0.7))            |      0.00694299 | 0.0 Mb      |         499 |         998 |                   0 |                  0 |  0.00305184  |
+--------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| Schedule|Annealing((3.5, 0.7)) |      0.00741506 | 0.015625 Mb |         499 |         998 |                   0 |                  0 | -6.09643e-06 |
+--------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

По результатам видно, что schedule версия достигла большей точности за то же число итераций. Значит применение schedule более чем осмыслено, поэтому в дальнейшем рассмотрении будет использоваться именно schedule версия.

Теперь рассмотрим метод имитации отжига на какой-нибудь функции, на которой легко сойтись, чтобы проверить эффектвность такого метода в тех ситуациях, когда классические методы показывают себя неплохо.

seed(1717701126228002400)

methods = [
    pow_step_schedule(5, 0.0001, 4, SimulatedAnnealing)(5, 0.3, 0.05),
    CoordinateDescent(learning_rate=3, aprox_dec=1e-30),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(100),
    AbsolutePrecisionCondition(1e-7, np.array([0]))
]

easy_function = LambdaOracul(lambda x: x**2)

result = Runner.run(methods, [easy_function], np.array([9.0]), conditions + modules, precision=1e-7,
                    **TABLE, **VISUALIZE, animate=True, animate_main=True, animation_main_full=True,
                    animator_main_step=0.01)

SEED: 1717701126228002400

LambdaOracul
+------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------+
| Method name                  |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount |      Result |
+==============================+=================+=============+=============+=============+=====================+====================+=============+
| Schedule|Annealing((5, 0.3)) |     0.001899    | 0.0 Mb      |          99 |         198 |                   0 |                  0 | 7.65804e-06 |
+------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------+
| CoordinateDescent(3)         |     0.000398636 | 0.015625 Mb |           3 |           7 |                   0 |                  0 | 0           |
+------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------+
| ScipyMethod(newton-cg)       |     0.00269103  | 0.09375 Mb  |           1 |           2 |                   2 |                  2 | 2.13163e-14 |
+------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Метод имитации отжига сходился медленее и менее точно чем остальные, что не удивительно зная суть данного алгоритма. На самом деле разница просто огромна, пока классическимЪ методами нужна всего пара итераций чтобы сойтись с невероятной точностью, методу имитации отжига нужны сотни итераций и это с учётом schedule, который тоже требовал первоночального подгона и который нехило увеличил точность по сравнению с обычным результатом. Исходя из вышесказанного можно сделать вывод, что использование метода имитации обжига на унимодальных функциях нецелесообразно.

Эволюционные методы оптимизации#

Недостаток метода имитации отжига виден на лицо, но можно ли его как-то исправить? Если каждая итерация этого метода такая легковесная, то может просто увеличить их количество? Логическим продолжением метода отжига можно назвать эволюционные алгоритмы. Основная суть таких алгоритмов - генерация потомков из родителей и использование лучших из представителей для следующего поколения.

(μ,λ)-ES и (μ+λ)-ES#

Сначала рассмотрим стандартные эволюционные методы: (μ,λ)-ES и (μ+λ)-ES. Алгоритм этих методов очень прост:

  1. Случайно относительно родителей генерируются потомки.

  2. Те представители, которые оказались лучше по характеристикам (в нашем случае, где значение оракула наименьшее) становятся родителями для следующего поколения.

  3. Повторяем

Отличие методов заключается лишь в том, что первый метод использует только потомков в следуюшем поколении, а второй и потомков и родителей.

В нашем случае они реализованы через один класс, имеющим такой алгоритм: На первом шагу случайно относительно стартовой точки генерируется необходимое количество родителей. На следующих шагах каждый ребенок генерируется относительно случайного родителя, среди всех потомков и какого-то числа лучших представителей родителей (данное число задано параметром, что и позволило через один алгоритм реализовать оба метода: в одном случае этот параметр равен 0, в другом числу родителей) выбираются лучшие, они станут родителями для следующего поколения. Точкой которая вовзращается на итерации - лучшая точка среди тех, что станут родителями на следующей итерации.

Методы параметризуются с помощью learning rate, максимального отклонения от родителя по каждой из осей, parents, количества родиетелей, childrens, количества детей, и reuse parents, числа лучших представителей родителей, которые могут использоваться в следующей итерации, знаение меньше 0 означает использование всех родителей.

Сравним эти методы на сложной функции из первого сравнения:

seed(1717703580539391700)

methods = [
    pow_step_schedule(3.5, 0.000001, 3, SimulatedAnnealing)(3.5, 0.7, 0.02),
    Evolution(50, 10, 50, 0),
    Evolution(50, 10, 50, 50),
    CoordinateDescent(learning_rate=3, aprox_dec=1e-30),
]

conditions = [
    StepCountCondition(500)
]

result = Runner.run(methods, [hard_function], np.array([9.0]), conditions + modules, precision=1e-7,
                    **TABLE, **VISUALIZE, animate=True, animate_main=True, animation_main_full=True,
                    animator_main_step=0.01)

SEED: 1717703580539391700

LambdaOracul
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| Method name                                            |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount |       Result |
+========================================================+=================+=============+=============+=============+=====================+====================+==============+
| Schedule|Annealing((3.5, 0.7))                         |      0.00735998 | 0.0 Mb      |         499 |         998 |                   0 |                  0 |  0.000610314 |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=0)  |      0.0851898  | 0.015625 Mb |         499 |       29950 |                   0 |                  0 |  0.461961    |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=50) |      0.0921309  | 0.0 Mb      |         499 |       34940 |                   0 |                  0 | -0.000101975 |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+
| CoordinateDescent(3)                                   |      0.00501513 | 0.03125 Mb  |         499 |         240 |                   0 |                  0 | 11.6239      |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+--------------+

Оба метода показали сходимость к глобальному минимуму, но (μ+λ)-ES достиг намного большей точности чем (μ,λ)-ES и даже обогнал метод имитации отжига, но кажется что метод может достигнуть большей точности если уменьшать learning rate с помощью schedule, но об этом позже. Превосходство метода, который не отбрасывает родителей очевиден - во-первых забывать о родителях это не по-христиански и бог будет всячески мешать, а во вторых возможна ситуация когда при генерации поколения были выбраны плохие родители или потомки просто сгенерировались неудачно, тогда использование только детей в следующем состоянии может ухудшить текущее положение, но такой подход может позволить вылезти из локальных минимумов. Также не стоит забывать, что итоговая точка - лучшая точка на последней итераци, поэтому возможна ситуация, когда какой-то из родителей уже находится возле глобального минимума, но ни один из потомков не сгенерировался возле него, тогда на следующей итерации об этой точке забудут и в зачёт пойдёт какой-то из менее успешных потомков.

Кажется, что главное преимущество эволюционных алгоритмов - множество ‘ветвей’ которые могут рассматриваться алгоритмом одновременно, на двумерном пространстве данное преимущество теряется поэтому рассмотрим функцию большей размерности.

В качестве функции для сравнения рассмотрим функцию леви из 1 лабораторной работы. При проведении эксперимента с этой функции все методы сошлись к локальным минимумам, которые достаточно далеко находились от глобального.

seed(1717709095419857600)

methods = [
    pow_step_schedule(70, 0.000001, 2, SimulatedAnnealing)(70, 0.7, 0.02),
    Evolution(50, 10, 50, 0),
    Evolution(50, 10, 50, -1),
    GradientDescent(learning_rate=10, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-2, [5, 5])
]


def levi_function(x, y):
    FACTOR = 5
    return (np.sin(3 / FACTOR * np.pi * x) ** 2
            + (x / FACTOR - 1) ** 2 * (1 + np.sin(3 * np.pi * y / FACTOR) ** 2)
            + (y / FACTOR - 1) ** 2 * (1 + np.sin(2 * np.pi * y / FACTOR) ** 2))


levi_oracul = LambdaOracul(levi_function)

point = np.array([1000.0, 1000.0])
result = Runner.run(methods, [levi_oracul], point, modules + conditions, precision=1e-7,
                    **TABLE, **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717709095419857600

/Users/mike/Documents/Methods-of-optimization/presentation/lab-4/../../src/v2/visualization/animation.py:76: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  figure_contour = plt.figure()
/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

LambdaOracul
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| Method name                                            |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                      |
+========================================================+=================+=============+=============+=============+=====================+====================+=============================+
| Schedule|Annealing((70, 0.7))                          |       0.0401847 | 0.015625 Mb |        1999 |        3998 |                   0 |                  0 | [6.60595672 3.42536351]     |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=0)  |       0.523847  | 0.0 Mb      |        1999 |      119950 |                   0 |                  0 | [7.02246432 8.94924154]     |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=10) |       0.579792  | 0.0 Mb      |        1999 |      139940 |                   0 |                  0 | [4.99035148 5.0193696 ]     |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| GradientDescent(10,GoldenRatioMethod,eps=1e-05)        |       1.41681   | 0.0 Mb      |        1999 |       61969 |                1999 |                  0 | [ -3.68941723 474.99862478] |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| ScipyMethod(newton-cg)                                 |       0.499323  | 0.0 Mb      |        1999 |         425 |                 425 |                400 | [287.23759857 999.99844021] |
+--------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+

Классические методы повторили свой результат - они сошлись к локальным минимумам, что ожидаемо, ведь данные методы не гарантируют сходимости на не унимодальных функциях. Стахастические методы показали себя намного лучше - все оказались возле глобального минимума, но (μ+λ)-ES показал наивысшую точность.

Сходимость методов легко объяснить, если рассмотреть общий вид функции на графиках. Если убрать выступы, то заметно сильное убываение функции по обеим осям в сторону глобального минимума. Стахастическим методам этого достаточно, чтобы с большим шансом двигаться именно в сторону глобального минимума и в итоге оказаться совсем близко к нему. При этом выступы почти никак не мешают методам в случае если learning rate больше ширины выступа, так как методы могут с легкостью перепрыгнуть через них.

Теперь хотелось бы рассмотреть эволюционные методы на унимодальных функциях, чтобы узнать, унаследовали ли он главный недостаток метода имитации обжига. В качестве функции для сравнения будет использовться функция Бута из 1 лабораторной работы.

seed(1717712609050922700)

methods = [
    pow_step_schedule(70, 0.0001, 2, SimulatedAnnealing)(70, 0.7, 0.02),
    Evolution(50, 10, 50, 0),
    Evolution(50, 10, 50, -1),
    CoordinateDescent(learning_rate=10, aprox_dec=1e-3),
    GradientDescent(learning_rate=100, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-5, [1, 3])
]

bute = LambdaOracul(lambda x, y: (x + 2 * y - 7) ** 2 + (2 * x + y - 5) ** 2)

result = Runner.run(methods, [bute], np.array([100, 100]), modules + conditions,
                    precision=1e-7, **TABLE, **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717712609050922700

LambdaOracul
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| Method name                                            |   ExecutionTime | RAMSize   |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                    |
+========================================================+=================+===========+=============+=============+=====================+====================+===========================+
| Schedule|Annealing((70, 0.7))                          |      0.0341501  | 0.0 Mb    |        1999 |        3998 |                   0 |                  0 | [0.99009294 2.99960765]   |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=0)  |      0.354077   | 0.0 Mb    |        1999 |      119950 |                   0 |                  0 | [-2.8398828   5.05916182] |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=10) |      0.372421   | 0.0 Mb    |        1999 |      139940 |                   0 |                  0 | [1.30833359 2.79108948]   |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| CoordinateDescent(10)                                  |      0.0216882  | 0.0 Mb    |        1999 |         105 |                   0 |                  0 | [1 2]                     |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| GradientDescent(100,GoldenRatioMethod,eps=1e-05)       |      0.00492692 | 0.0 Mb    |           5 |         180 |                   5 |                  0 | [1.00000766 3.00000103]   |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| ScipyMethod(newton-cg)                                 |      0.00622797 | 0.0 Mb    |           1 |           5 |                   5 |                  5 | [1. 3.]                   |
+--------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Градиентный спуск и Newton-CG ожидаемо сошлись, но к удивлению координатный спуск сошелся в точку [1, 2], хотя минимум был в [1, 3]. Стахастические методы тоже сошлись к глобальному минимуму, причем метод имитации отжига сошелся достаточно близко, в то время как эволюционные методы показали очень низкую точность. Данный недостаток можно решить с помощью schedule, но существует более практичное решение - это логическое продолжение эволюционных методов.

CMA-ES#

Cледующим методом после (μ,λ)-ES и (μ+λ)-ES является CMA-ES. Это метод адаптации матриц ковариации. Суть этого метода почти полностью совпадает с вышеописанными эволюционными методами за исключением того, как потомки генерируются из родителей. Для этого на основе родителей находится мат ожидание и матрица ковариации, с помощью которых строится нормальное распределение, которое с наибольшим шансом могло сгенерировать данную конфигурацию точек. С помощью найденного распределения случайным образом генерируются потомки. В остальном метод совпадает с другими эволюционными методами. Главный неочевидный плюс CMA-ES заключается в распределении, которое получается на каждой итерации. Если метод уже находится возле глобального минимума и все родители находятся недалеко друг от друга, то генерируется распределение с небольшим разбросом, то есть потомки будут генерироваться в небольшой области. Другими словами, по мере приближения к глобальному минимуму разброс будет уменьшаться и проблема с необходимостью уменьшения learning rate не возникнет.

Метод параметризуeтся с помощью learning rate, максимального отклонения по каждой из осей от стартовой точки для первого поколения, parents, количества родиетелей, childrens, количества детей, и reuse parents, числа лучших представителей родителей, которые могут использоваться в следующей итерации, знаение меньше 0 означает использование всех родителей.

Сначала проверим работоспомобность на сложной функции, например на функции Экли.

seed(1717716266320449600)

methods = [
    pow_step_schedule(70, 0.00001, 2, SimulatedAnnealing)(70, 0.3, 0.005),
    Evolution(50, 10, 50, -1),
    EvolutionCMA(200, 10, 100, -1),
    GradientDescent(learning_rate=50, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-5, [0, 0])
]

ekli = LambdaOracul(
    lambda x, y: -20 * exp(-0.2 * sqrt((x ** 2 + y ** 2) / 2)) - exp((cos(2 * pi * x) + cos(2 * pi * y)) / 2) + e + 20)

result = Runner.run(methods, [ekli], np.array([100, 100]), modules + conditions, precision=1e-7, **TABLE, **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717716266320449600

LambdaOracul
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+
| Method name                                                 |   ExecutionTime | RAMSize   |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                          |
+=============================================================+=================+===========+=============+=============+=====================+====================+=================================+
| Schedule|Annealing((70, 0.3))                               |      0.0346618  | 0.0 Mb    |        1999 |        3998 |                   0 |                  0 | [126.99904171  -0.99327477]     |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=10)      |      0.394989   | 0.0 Mb    |        1999 |      139940 |                   0 |                  0 | [-0.1434135  -0.01287891]       |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+
| EvolutionCMA(200,parents=10,childrens=100,reuse_parents=10) |      0.00402498 | 0.0 Mb    |          11 |        1340 |                   0 |                  0 | [6.30016662e-06 6.51497808e-06] |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+
| GradientDescent(50,GoldenRatioMethod,eps=1e-05)             |      1.22463    | 0.125 Mb  |        1999 |       69965 |                1999 |                  0 | [18.99916573 18.99916573]       |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+
| ScipyMethod(newton-cg)                                      |      0.0252311  | 0.0 Mb    |        1999 |           2 |                   2 |                  1 | [100. 100.]                     |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Классические методы ожидаемо не сошлись. Стахастические повели себя по-разному. Метод имитации отжига не сошелся к глобальному минимуму, в отличие от функции леви на данной функции нет заметного спада в сторону глобального минимума, поэтому локальные минимумы почти равносильны для данного метода и схождение к глобальному минимуму может быть не более чем удачей. Эволюционные методы же сошлись, как раз из-за большего количества локальных минимумов, которые может рассматривать метод. (μ+λ)-ES сошелся с низкой точностью и довольно медлено из-за описанных ранее проблем, в то время как CMA-ES сошелся невероятно быстро и невероятно точно. Трудно теоретически объяснить почему функция показала такой отличный результат.

Сравним эффективность данного метода на унимодальной функции, например на функции Химмельблау.

seed(1717717613926377900)
methods = [
    pow_step_schedule(70, 0.00001, 2, SimulatedAnnealing)(70, 0.3, 0.005),
    Evolution(50, 10, 50, -1),
    EvolutionCMA(200, 10, 100, -1),
    GradientDescent(learning_rate=50, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]

himmelblau = LambdaOracul(lambda x, y: (x ** 2 + y - 11) ** 2 + (x + y ** 2 - 7) ** 2)

conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-2, [3, 2]),
    AbsolutePrecisionCondition(1e-2, [-2.805118, 3.131312]),
    AbsolutePrecisionCondition(1e-2, [-3.779310, -3.283186]),
    AbsolutePrecisionCondition(1e-2, [3.584428, -1.8481126]),
    PrecisionCondition(1e-4)]

result = Runner.run(methods, [himmelblau], np.array([100, 100]), modules + conditions, precision=1e-7, **TABLE,
                    **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717717613926377900

LambdaOracul
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| Method name                                                 |   ExecutionTime | RAMSize   |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                    |
+=============================================================+=================+===========+=============+=============+=====================+====================+===========================+
| Schedule|Annealing((70, 0.3))                               |      0.0506401  | 0.0 Mb    |        1999 |        3998 |                   0 |                  0 | [-2.90318959  3.18086686] |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=10)      |      0.387323   | 0.0 Mb    |        1999 |      139940 |                   0 |                  0 | [2.96313865 2.0214291 ]   |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| EvolutionCMA(200,parents=10,childrens=100,reuse_parents=10) |      0.00494909 | 0.0 Mb    |          17 |        2060 |                   0 |                  0 | [3.00365622 1.99654897]   |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| GradientDescent(50,GoldenRatioMethod,eps=1e-05)             |      0.00716996 | 0.0 Mb    |           8 |         280 |                   8 |                  0 | [-3.77519692 -3.27807661] |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+
| ScipyMethod(newton-cg)                                      |      0.016875   | 0.0 Mb    |           1 |          17 |                  17 |                 17 | [3. 2.]                   |
+-------------------------------------------------------------+-----------------+-----------+-------------+-------------+---------------------+--------------------+---------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Здесь CMA-ES показал себя отлично, всего на пару итераций медленее сойдясь чем метод градиентного спуска. Показав схожие результаты по времени исполнения и требуемой оперативной памяти, он сделал кратно больше вызовов. То есть в случае, если вызов оракула может быть тяжеловестной операцией, как например во многих задачах машиного обучения, то CMA-ES покажет результат намного хуже, но не смотря на это CMA-ES решает намного большее количество задач чем могут решить классические алгоритмы.

Общее сравнение#

Определившись с плюсами и минусами описанных методов, хотелось бы провести ещё несколько тестов для подтверждения выдвинутых утверждений, объясняющих поведение методов на тех или иных функциях.

Начнём с функции Розенброка, в предыдущих лабараторных работах классические методы показали себя относительной плохо, исключением является Newton-CG, который показал отличные результаты.

seed(1717718715696515100)

methods = [
    pow_step_schedule(20, 0.00001, 2, SimulatedAnnealing)(70, 0.3, 0.005),
    Evolution(10, 10, 50, -1),
    EvolutionCMA(100, 10, 100, -1),
    GradientDescent(learning_rate=10, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-5, [1, 1])
]
rosenbrok = LambdaOracul(lambda x, y: (1 - np.float64(x)) ** 2 + 100 * (np.float64(y) - np.float64(x) ** 2) ** 2)

result = Runner.run(methods, [rosenbrok], np.array([50, -50]), modules + conditions,
                    precision=1e-7, **TABLE, **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717718715696515100

LambdaOracul
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+
| Method name                                                 |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                  |
+=============================================================+=================+=============+=============+=============+=====================+====================+=========================+
| Schedule|Annealing((70, 0.3))                               |      0.0344238  | 0.0 Mb      |        1999 |        3998 |                   0 |                  0 | [1.21141131 1.46872273] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+
| Evolution(10,parents=10,childrens=50,reuse_parents=10)      |      0.392035   | 0.0 Mb      |        1999 |      139940 |                   0 |                  0 | [1.04019468 1.073527  ] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+
| EvolutionCMA(100,parents=10,childrens=100,reuse_parents=10) |      0.00485802 | 0.0 Mb      |          15 |        1820 |                   0 |                  0 | [0.99999514 0.99999566] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+
| GradientDescent(10,GoldenRatioMethod,eps=1e-05)             |      1.1518     | 0.046875 Mb |        1999 |       61969 |                1999 |                  0 | [0.9074093  0.82321555] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+
| ScipyMethod(newton-cg)                                      |      0.0464041  | 0.0 Mb      |           1 |          53 |                  53 |                 42 | [1. 1.]                 |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Результат эксперимента из предыдущих лабораторных работ повторился - Newton-CG сошелся, градиентный спуск не сошелся, но оказался довольно близко. Стахастические методы показали результат похожий на предыдущие тесты, метод имитации отжига был близок, но не точен и медленен, (μ+λ)-ES был быстрее и точнее, но все равно недостаточно, CMA-ES всего за немного больше десятка операций сошелся с удивительной точностью.

Проверим CMA-ES на функции с которой никто не справлялся - Изома.

seed(1717719845694252000)

methods = [
    pow_step_schedule(70, 0.00001, 2, SimulatedAnnealing)(70, 0.3, 0.005),
    Evolution(50, 10, 50, -1),
    EvolutionCMA(200, 10, 100, -1),
    GradientDescent(learning_rate=50, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-5, [pi, pi])
]

isom = LambdaOracul(lambda x, y: -np.cos(x) * np.cos(y) * np.exp(-((x - np.pi) ** 2 + (y - np.pi) ** 2)))

result = Runner.run(methods, [isom], np.array([50, 50]), modules + conditions,
                    precision=1e-7, **TABLE, **VISUALIZE, **FULL_VISUALIZE)

SEED: 1717719845694252000

LambdaOracul
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| Method name                                                 |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                      |
+=============================================================+=================+=============+=============+=============+=====================+====================+=============================+
| Schedule|Annealing((70, 0.3))                               |      0.0381141  | 0.0 Mb      |        1999 |        3998 |                   0 |                  0 | [161.30771744 294.26458551] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| Evolution(50,parents=10,childrens=50,reuse_parents=10)      |      0.512727   | 0.0 Mb      |        1999 |      139940 |                   0 |                  0 | [3.27860234 3.08610067]     |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| EvolutionCMA(200,parents=10,childrens=100,reuse_parents=10) |      0.00544906 | 0.015625 Mb |          13 |        1580 |                   0 |                  0 | [3.14158831 3.14160154]     |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| GradientDescent(50,GoldenRatioMethod,eps=1e-05)             |      1.37478    | 0.09375 Mb  |        1999 |       69965 |                1999 |                  0 | [50. 50.]                   |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+
| ScipyMethod(newton-cg)                                      |      0.0234089  | 0.0 Mb      |        1999 |           1 |                   1 |                  1 | [50. 50.]                   |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Здесь сошлись только (μ+λ)-ES и CMA-ES, но второй опять сделал это быстрее и намного точнее.

Пришло время для финального босса - авторской функции, которую было решено назвать функцией Скорнякова. Это функция с множеством локальных минимумов, и даже недалеко от глобального минимума есть множество локальных минимумов, которые должны стать ловушкой для стахастических методов.

seed()

methods = [
    pow_step_schedule(50, 0.00001, 2, SimulatedAnnealing)(70, 0.3, 0.005),
    Evolution(10, 10, 50, -1),
    EvolutionCMA(200, 10, 100, -1),
    GradientDescent(learning_rate=50, aprox_dec=1e-5),
    ScipyMethod("Newton-CG")
]
conditions = [
    StepCountCondition(2000),
    AbsolutePrecisionCondition(1e-5, [0, 0])
]

skornyakov_function = LambdaOracul(lambda x, y: sqrt((x + 4 * sin(x)) ** 2 + (y + 2 * sin(y)) ** 2))

result = Runner.run(methods, [skornyakov_function], np.array([50, -50]), modules + conditions, precision=1e-7, **TABLE, **VISUALIZE,**FULL_VISUALIZE)

SEED: 1718012026564032000

LambdaOracul
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+
| Method name                                                 |   ExecutionTime | RAMSize     |   StepCount |   CallCount |   GradientCallCount |   HessianCallCount | Result                            |
+=============================================================+=================+=============+=============+=============+=====================+====================+===================================+
| Schedule|Annealing((70, 0.3))                               |      0.0338032  | 0.0 Mb      |        1999 |        3998 |                   0 |                  0 | [0.00796779 0.02188699]           |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+
| Evolution(10,parents=10,childrens=50,reuse_parents=10)      |      0.373072   | 0.0 Mb      |        1999 |      139940 |                   0 |                  0 | [-0.00431426 -0.03712353]         |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+
| EvolutionCMA(200,parents=10,childrens=100,reuse_parents=10) |      0.00432324 | 0.0 Mb      |          15 |        1820 |                   0 |                  0 | [ 4.50640336e-07 -2.33714497e-06] |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+
| GradientDescent(50,GoldenRatioMethod,eps=1e-05)             |      1.20276    | 0.078125 Mb |        1999 |       69965 |                1999 |                  0 | [-10.737195   -10.45537135]       |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+
| ScipyMethod(newton-cg)                                      |      0.03022    | 0.0 Mb      |        1999 |           9 |                   9 |                  7 | [ 42.15882057 -48.17108736]       |
+-------------------------------------------------------------+-----------------+-------------+-------------+-------------+---------------------+--------------------+-----------------------------------+

/Users/mike/Documents/Methods-of-optimization/venv/lib/python3.9/site-packages/matplotlib/animation.py:892: UserWarning: Animation was deleted without rendering anything. This is most likely not intended. To prevent deletion, assign the Animation to a variable, e.g. `anim`, that exists until you output the Animation using `plt.show()` or `anim.save()`.
  warnings.warn(

Результат не удивил, классические функции разумеется не сошлись. Метод имитации отжига попал в один из локальных минимумов из которого не смог вылезти, поэтому он не сошелся. Эволюционные методы сошлись и опять CMA-ES сделал это очень быстро и с высокой точностью.

Contents